Sys.setenv(RETICULATE_PYTHON="/home/fdeckert/bin/miniconda3/envs/p.3.8.12-FD20200109SPLENO/bin/python3")
library_load <- suppressMessages(
list(
# Seurat
library(Seurat),
library(SeuratWrappers),
library(monocle3),
# RaceID
require(Matrix),
require(RaceID),
# SingleR
library(SingleR),
library(SingleCellExperiment),
# TradeSeq
library(tradeSeq),
library(clusterExperiment),
# GSEA
library(msigdbr),
library(fgsea),
# Data
library(tidyverse),
library(openxlsx),
# Plotting
library(ComplexHeatmap),
library(circlize),
library(viridis),
library(ggplotify),
# Python
library(reticulate)
)
)
random_seed <- 42
set.seed(random_seed)
options(warn=-1)
ht_opt$message=FALSE # ComplexHeatmap
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
# Source files
source("plotting_global.R")
source("bin/cell_type.R")
source("bin/seurat_qc.R")
source("bin/seurat_dea.R")
source("bin/cell_type.R")
source("bin/raceid.R")
source("bin/gene_modules.R")
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()
so_file <- "data/object/sct/so_sct_int_hvg8000.rds"
so_pp_file <- "data/object/pp.rds"
h5ad_pp_file <- "data/object/pp.h5ad"
so <- readRDS(so_file)
DefaultAssay(so) <- "integrated"
# Cluster all cells
so <- FindNeighbors(so, dims=1:10, k.param=20, verbose=FALSE)
so <- FindClusters(so, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
# Umap dimensional reduction
so <- RunUMAP(so, dims=1:15, n.neighbors=100, min.dist=1, spread=1, verbose=FALSE)
options(repr.plot.width=30, repr.plot.height=10)
dplot_1 <- dplot(so, reduction="umap", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so, reduction="umap", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_3 <- dplot(so, reduction="umap", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
fplot_1 <- fplot(so, reduction="umap", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
fplot_2 <- fplot(so, reduction="umap", features="pMt_RNA") + ggtitle("Percentage Mt") + scale_color_viridis(option="G")
fplot_3 <- fplot(so, reduction="umap", features="Ccr2") + ggtitle("Ccr2") + scale_color_viridis(option="G")
fplot_4 <- fplot(so, reduction="umap", features="Ccr7") + ggtitle("Ccr7") + scale_color_viridis(option="G")
fplot_5 <- fplot(so, reduction="umap", features="Ly6c2") + ggtitle("Ly6c2") + scale_color_viridis(option="G")
dplot_1 + dplot_2 + dplot_3 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(nrow=2, ncol=6)
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
so@meta.data %>% dplyr::group_by(seurat_clusters, treatment, sample_rep) %>%
dplyr::summarise(n=n()) %>% data.frame() %>%
tidyr::spread(seurat_clusters, n) %>%
kableExtra::kable("html") %>% as.character() %>% IRdisplay::display_html()
`summarise()` has grouped output by 'seurat_clusters', 'treatment'. You can override using the `.groups` argument.
| treatment | sample_rep | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CpG | Rep1 | 705 | 590 | 154 | 476 | 412 | 167 | 383 | 337 | 258 | 254 | 157 | 154 | 102 | 56 | 133 | 117 | 100 | 70 | 101 | 74 | 13 | NA |
| CpG | Rep2 | 930 | 654 | 160 | 648 | 108 | 348 | 106 | 574 | 226 | 450 | 343 | 314 | 296 | 150 | 293 | 98 | 66 | 161 | 38 | 20 | 43 | 36 |
| NaCl | Rep1 | 170 | 161 | 327 | 72 | 372 | 106 | 437 | 23 | 321 | 49 | 121 | 77 | 59 | 133 | 70 | 65 | 121 | 42 | 106 | 109 | 15 | NA |
| NaCl | Rep2 | 257 | 251 | 685 | 90 | 224 | 477 | 91 | 52 | 138 | 167 | 252 | 303 | 314 | 426 | 117 | 73 | 60 | 61 | 31 | 33 | 25 | 4 |
# Rmove cluster wich have no evidence in replicates
so <- subset(so, subset=seurat_clusters!="21")
DefaultAssay(so) <- "RNA"
raceid <- raceid_pp(so=so, suffix="", compute=FALSE)
so <- raceid_to_seurat(so, raceid)
options(repr.plot.width=30, repr.plot.height=5)
dplot_1 <- dplot(so, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
dplot_5 <- dplot(so, reduction="umap_varid", group_by="label_fine_haemosphere", alpha=0.5) + scale_color_manual(values=color$label_fine_haemosphere)
fplot_1 <- fplot(so, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + plot_layout(nrow=1, ncol=6)
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
so@meta.data %>% dplyr::group_by(varid_clusters, treatment, sample_rep) %>%
dplyr::summarise(n=n()) %>% data.frame() %>%
tidyr::spread(varid_clusters, n) %>%
kableExtra::kable("html") %>% as.character() %>% IRdisplay::display_html()
`summarise()` has grouped output by 'varid_clusters', 'treatment'. You can override using the `.groups` argument.
| treatment | sample_rep | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| CpG | Rep1 | 904 | 270 | 682 | 155 | 565 | 519 | 464 | 326 | 157 | 195 | 186 | 89 | 67 | 45 | 59 | 41 | 31 | 42 | 12 | 3 | 1 |
| CpG | Rep2 | 1294 | 653 | 984 | 207 | 266 | 262 | 664 | 447 | 340 | 183 | 277 | 64 | 83 | 94 | 29 | 88 | 62 | 10 | 13 | 5 | 1 |
| NaCl | Rep1 | 153 | 204 | 116 | 404 | 551 | 458 | 155 | 49 | 118 | 305 | 73 | 125 | 37 | 41 | 65 | 10 | 26 | 27 | 11 | 19 | 9 |
| NaCl | Rep2 | 445 | 935 | 235 | 893 | 171 | 306 | 209 | 107 | 246 | 101 | 120 | 61 | 47 | 49 | 28 | 41 | 37 | 41 | 26 | 16 | 13 |
seurat_clusters_prog <- c("17", "15", "2", "13", "5", "12", "11", "9", "7", "3", "0", "1")
so_prog <- subset(so, subset=seurat_clusters %in% seurat_clusters_prog)
so_prog <- AddModuleScore(so_prog, list(genes_hsc), assay="RNA", name="msHSC")
so_prog <- AddModuleScore(so_prog, list(genes_ly), assay="RNA", name="msLy")
so_prog <- AddModuleScore(so_prog, list(genes_meg), assay="RNA", name="msMeg")
so_prog <- AddModuleScore(so_prog, list(genes_ery), assay="RNA", name="msEry")
so_prog <- AddModuleScore(so_prog, list(genes_mo), assay="RNA", name="msMo")
so_prog <- AddModuleScore(so_prog, list(genes_neu), assay="RNA", name="msNeu")
so_prog <- AddModuleScore(so_prog, list(genes_eo), assay="RNA", name="msEo")
so_prog <- AddModuleScore(so_prog, list(genes_baso), assay="RNA", name="msBaso")
so_prog <- AddModuleScore(so_prog, list(genes_mast), assay="RNA", name="msMast")
raceid_prog <- raceid_pp(so=so_prog, suffix="_prog", compute=FALSE)
so_prog <- raceid_to_seurat(so_prog, raceid_prog)
options(repr.plot.width=30, repr.plot.height=5)
dplot_1 <- dplot(so_prog, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so_prog, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so_prog, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so_prog, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
dplot_5 <- dplot(so_prog, reduction="umap_varid", group_by="label_fine_haemosphere", alpha=0.5) + scale_color_manual(values=color$label_fine_haemosphere)
fplot_1 <- fplot(so_prog, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + plot_layout(nrow=1, ncol=6)
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
options(repr.plot.width=7.5, repr.plot.height=10)
# Get expression matrix
mat <- AverageExpression(so_prog, assay="RNA", slot="data", features=genes_marker$genes, group.by="seurat_clusters")[["RNA"]]
mat <- mat[, seurat_clusters_prog]
mat <- t(scale(t(mat)))
breaks <- seq(-max(abs(mat)), +max(abs(mat)), by=0.1)
hm_1 <- ComplexHeatmap::Heatmap(
matrix=mat,
name="z-score",
column_title="Progenitor marker expression",
col=colorRamp2(breaks, mako(length(breaks))),
column_title_gp=gpar(fontsize=10, fontface="bold"),
column_names_gp =grid::gpar(fontsize=8),
row_title_gp=gpar(fontsize=8, fontface="bold"),
row_names_gp=grid::gpar(fontsize=8),
row_split=genes_marker$cell_type,
cluster_rows=TRUE,
cluster_columns=FALSE,
show_row_names=TRUE,
show_column_names=TRUE,
row_dend_width=unit(1.5, "cm"),
width=ncol(mat)*unit(3, "mm"),
height=nrow(mat)*unit(6, "mm"),
rect_gp=gpar(col="white", lwd=2),
heatmap_legend_param=list(title_gp=gpar(fontsize=8, fontface="bold"), labels_gp=gpar(fontsize=8))
) %>% as.ggplot()
hm_1
seurat_clusters_eb <- c("13", "5", "12", "11", "9", "7", "3", "0", "1")
so_eb <- subset(so, subset=seurat_clusters %in% seurat_clusters_eb)
cell_type_eb <- data.frame(
ident=c(
"13",
"5",
"12",
"11",
"9",
"7",
"3",
"0",
"1"
),
cell_type_fine_detail=c(
"ProEB (1)",
"ProEB (2)",
"ProEB (3)",
"ProEB (4)",
"EB (1)",
"EB (2)",
"EB (3)",
"EB (4)",
"EB (5)"
),
cell_type_fine=c(
"ProEB (1)",
"ProEB (2)",
"ProEB (3)",
"ProEB (4)",
"EB (1)",
"EB (2)",
"EB (3)",
"EB (4)",
"EB (5)"
),
cell_type_main=c(
"ProEB",
"ProEB",
"ProEB",
"ProEB",
"EB",
"EB",
"EB",
"EB",
"EB"
)
)
cell_type_eb <- dplyr::left_join(dplyr::select(so_eb@meta.data, seurat_clusters, cell_id), cell_type_eb, by=c("seurat_clusters"="ident")) %>%
column_to_rownames("cell_id") %>%
dplyr::mutate(varid_clusters=paste0("varid_clusters_eb_NA")) %>%
dplyr::mutate(ident=paste0("seurat_clusters_eb_", seurat_clusters))
seurat_clusters_mpp <- c("17", "15", "2")
so_mpp <- subset(so_prog, subset=seurat_clusters %in% seurat_clusters_mpp)
raceid_mpp <- raceid_pp(so=so_mpp, suffix="_mpp", compute=FALSE)
so_mpp <- raceid_to_seurat(so_mpp, raceid_mpp)
options(repr.plot.width=30, repr.plot.height=15)
dplot_1 <- dplot(so_mpp, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so_mpp, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so_mpp, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so_mpp, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
fplot1 <- fplot(so_mpp, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
fplot2 <- fplot(so_mpp, reduction="umap_varid", features="pMt_RNA") + ggtitle("Percentage Mt") + scale_color_viridis(option="G")
fplot3 <- fplot(so_mpp, reduction="umap_varid", features="msHSC1") + ggtitle("msHSC1") + scale_color_viridis(option="G")
fplot4 <- fplot(so_mpp, reduction="umap_varid", features="msLy1") + ggtitle("msLy1") + scale_color_viridis(option="G")
fplot5 <- fplot(so_mpp, reduction="umap_varid", features="msMeg1") + ggtitle("msMeg1") + scale_color_viridis(option="G")
fplot6 <- fplot(so_mpp, reduction="umap_varid", features="msEry1") + ggtitle("msEry1") + scale_color_viridis(option="G")
fplot7 <- fplot(so_mpp, reduction="umap_varid", features="msMo1") + ggtitle("msMo1") + scale_color_viridis(option="G")
fplot8 <- fplot(so_mpp, reduction="umap_varid", features="msNeu1") + ggtitle("msNeu1") + scale_color_viridis(option="G")
fplot9 <- fplot(so_mpp, reduction="umap_varid", features="msEo1") + ggtitle("msEo1") + scale_color_viridis(option="G")
fplot10 <- fplot(so_mpp, reduction="umap_varid", features="msBaso1") + ggtitle("msBaso1") + scale_color_viridis(option="G")
fplot11 <- fplot(so_mpp, reduction="umap_varid", features="msMast1") + ggtitle("msMast1") + scale_color_viridis(option="G")
dplot_1 + dplot_2 + dplot_3 + dplot_4 + fplot1 + fplot2 + fplot3 + fplot4 + fplot5 + fplot6 + fplot7 + fplot8 + fplot9 + fplot10 + fplot11 + plot_layout(ncol=6)
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
options(repr.plot.width=7.5, repr.plot.height=10)
source("bin/gene_modules.R")
mat <- AverageExpression(so_mpp, assay="RNA", slot="data", features=genes_marker$genes, group.by="varid_clusters")[["RNA"]]
mat <- t(scale(t(mat)))
mat <- na.omit(mat)
cell_type_mpp <- genes_marker$cell_type[genes_marker$genes %in% rownames(mat)]
breaks <- seq(-max(abs(mat)), +max(abs(mat)), by=0.1)
hm_1 <- ComplexHeatmap::Heatmap(
matrix=mat,
name="z-score",
column_title="Progenitor marker expression",
col=colorRamp2(breaks, mako(length(breaks))),
column_title_gp=gpar(fontsize=10, fontface="bold"),
column_names_gp =grid::gpar(fontsize=8),
row_title_gp=gpar(fontsize=8, fontface="bold"),
row_names_gp=grid::gpar(fontsize=8),
row_split=cell_type_mpp,
cluster_rows=TRUE,
cluster_columns=FALSE,
show_row_names=TRUE,
show_column_names=TRUE,
row_dend_width=unit(1.5, "cm"),
width=ncol(mat)*unit(3, "mm"),
height=nrow(mat)*unit(6, "mm"),
rect_gp=gpar(col="white", lwd=2),
heatmap_legend_param=list(title_gp=gpar(fontsize=8, fontface="bold"), labels_gp=gpar(fontsize=8))
) %>% as.ggplot()
hm_1
cell_type_mpp <- data.frame(
ident=c(
1,
2,
3,
4,
5,
6,
7,
8,
9,
10
),
cell_type_fine_detail=c(
"MEP (2)",
"MEP (4)",
"MEP (3)",
"GMP",
"MegP",
"BasoP",
"MLP",
"MDP",
"MEP (1)",
"MastP"
),
cell_type_fine=c(
"MEP (2)",
"MEP (4)",
"MEP (3)",
"GMP",
"MegP",
"BasoP",
"MLP",
"MDP",
"MEP (1)",
"MastP"
),
cell_type_main=c(
"MEP",
"MEP",
"MEP",
"GMP",
"MegP",
"GMP",
"MLP",
"MDP",
"MEP",
"GMP"
)
)
cell_type_mpp <- dplyr::left_join(dplyr::select(so_mpp@meta.data, seurat_clusters, varid_clusters, cell_id), cell_type_mpp, by=c("varid_clusters"="ident")) %>%
column_to_rownames("cell_id") %>%
dplyr::mutate(varid_clusters=paste0("varid_clusters_mpp_", varid_clusters)) %>%
dplyr::mutate(ident=varid_clusters)
seurat_clusters_m <- c("10", "6", "18", "20", "19", "4", "8", "14", "16")
so_m <- subset(so, subset=seurat_clusters %in% seurat_clusters_m)
raceid_m <- raceid_pp(so=so_m, suffix="_m", compute=FALSE)
so_m <- raceid_to_seurat(so_m, raceid_m)
options(repr.plot.width=30, repr.plot.height=7.5)
dplot_1 <- dplot(so_m, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so_m, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so_m, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so_m, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
fplot_1 <- fplot(so_m, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
fplot_2 <- fplot(so_m, reduction="umap_varid", features="pMt_RNA") + ggtitle("Percentage Mt") + scale_color_viridis(option="G")
dplot_1 + dplot_2 + dplot_3 + dplot_4 + fplot_1 + fplot_2 + plot_layout(nrow=1, ncol=6)
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
options(repr.plot.width=30, repr.plot.height=7.5)
fplot_1 <- fplot(so_m, reduction="umap_varid", features="Ccr7") + ggtitle("Ccr7 (TCZ)") + scale_color_viridis(option="G")
fplot_2 <- fplot(so_m, reduction="umap_varid", features="Cxcr5") + ggtitle("Cxcr5 (BCZ)") + scale_color_viridis(option="G")
fplot_3 <- fplot(so_m, reduction="umap_varid", features="Cxcr4") + ggtitle("Cxcr4 (RP)") + scale_color_viridis(option="G")
fplot_4 <- fplot(so_m, reduction="umap_varid", features="Ccr2") + ggtitle("Ccr2 (Blood)") + scale_color_viridis(option="G")
fplot_1 + fplot_2 + fplot_3 + fplot_4 + plot_spacer() + plot_spacer() + plot_layout(nrow=1, ncol=6)
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale. Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
so_m$varid_clusters <- paste0("varid_clusters_m_", so_m$varid_clusters)
features_m <- c(
"Adgre1", # F4/20
"Itgax", # Cd11c
"Itgam", # Cd11b
"Ly6c2", # Ly6c
"Csf1r", # Cd115
"Cd8a", # Cd8
"Cd4", # Cd4
"Napsa", # Napsa
"Lsp1" # Lsp1
)
options(repr.plot.width=2*(2.5+length(features_m)/2.5), repr.plot.height=5)
dp_1 <- dp_feature(so_m, features_m, group_by="varid_clusters", scale=FALSE, title="Myeloid marker")
dp_2 <- dp_feature(so_m, features_m, group_by="varid_clusters", scale=TRUE, title="Myeloid marker")
dp_1 + dp_2
`summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. `summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.
options(repr.plot.width=30, repr.plot.height=15)
varid_clusters_m <- c("varid_clusters_m_10", "varid_clusters_m_6", "varid_clusters_m_4", "varid_clusters_m_5", "varid_clusters_m_2", "varid_clusters_m_3", "varid_clusters_m_8", "varid_clusters_m_9", "varid_clusters_m_1", "varid_clusters_m_11", "varid_clusters_m_7")
dea_m <- dea_seurat(so_m, ident="varid_clusters", file="result/dea/dea_m", conserved=FALSE, only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
dea_m_conserved <- dea_seurat(so_m, ident="varid_clusters", file="result/dea/dea_m_conserved", conserved=TRUE, grouping_var="treatment", only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
hm_1 <- hm_dea(dea_m, so_m, column_name="varid_clusters", column_order=varid_clusters_m, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL)
hm_2 <- hm_dea(dea_m_conserved, so_m, column_name="varid_clusters", column_order=varid_clusters_m, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL, conserved=TRUE)
hm_1 + hm_2
features_mac <- c(
"Spic", # Spic | RPM marker
"Vcam1", # Vcam1 | RPM marker
"Slc40a1", # Slc40a1 | RPM marker | Iron transport
"Siglec1", # Cd169 | MMM marker. Binds vaious cell type.
"Nr1h3", # Lxr alpha | MMMM marker. A transcription factor activated by oxysterols
"Marco", # Marco | MZM marker
"Cd209b", # SIGNR1 | MZM marker
"Gas6", # Gas6 | TBM
"Mfge8"
)
split_mac <- c(
"RPM",
"RPM",
"RPM",
"MMM",
"MMM",
"MZM",
"MZM",
"TBM",
"TBM"
)
options(repr.plot.width=2*(2.5+length(features_mac)/2.5), repr.plot.height=5)
dp_1 <- dp_feature(so_m, features_mac, split_mac, group_by="varid_clusters", scale=FALSE, title="Macrophage marker")
dp_2 <- dp_feature(so_m, features_mac, split_mac, group_by="varid_clusters", scale=TRUE, title="Macrophage marker")
dp_1 + dp_2
`summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. `summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.
options(repr.plot.width=10, repr.plot.height=5)
varid_clusters_mac <- c("varid_clusters_m_10", "varid_clusters_m_6")
dea_mac <- dea_seurat(subset(so_m, subset=varid_clusters %in% varid_clusters_mac), ident="varid_clusters", file="result/dea/dea_mac", conserved=FALSE, only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
dea_mac_conserved <- dea_seurat(subset(so_m, subset=varid_clusters %in% varid_clusters_mac), ident="varid_clusters", file="result/dea/dea_mac_conserved", conserved=TRUE, grouping_var="treatment", only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
hm_1 <- hm_dea(dea_mac, subset(so_m, subset=varid_clusters %in% varid_clusters_mac), column_name="varid_clusters", column_order=varid_clusters_m, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL)
hm_2 <- hm_dea(dea_mac_conserved, subset(so_m, subset=varid_clusters %in% varid_clusters_mac), column_name="varid_clusters", column_order=varid_clusters_m, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL, conserved=TRUE)
hm_1 + hm_2
features_dc <- c(
"Zbtb46", # Zbtb46 | Common DC progenitor (CDP) | TF
"Xcr1", # Xcr1 | cDC1 of which most express CD8aa
"Ly75", # Dec205 (Lectin receptor) | cDC1 subset in the WP | After immunization cDC1 colocalize with and preferentially activate Cd8+ T cells in the WP
"Itgae", "Cd207", # Cd103 Langerin | cDC1 in the MZ and RP at stready-state
"Sirpa", "Gpr183", "Mbtps1", # Sirpa Ebi2 (recognizes oxysterol ligands) S1P | cDC2 in the bridging channel (BC) at steady-state condition and also express Cd11b. Ebi2 and S1p maintain BC localization.
"Esam", "Notch2", "Rbpj", "Clec4a4", "Irf4", "Klf4", # Esam (endothelial cell adhesion moleculoe) Notch2 Rbpj Dcir2 | cDC2 Esam hi subset | Notch2 and Rbpj signaling is required for Esam expression. The subset also express Cd11b, Cd4 and Dcir2
"Cx3cr1", "Il12a", "Il12b", "Tnf", # Cx3cr1 Il12 Il12 Infalpha | cDC2 Esam low subset | Independent of Notch2 and Irf4 signaling and producing innflammatory cytokines such as TNFa and Il12. Also positve for Cd11b, Dcir2 but less cd4 or double negative.
"Cd4",
"Cd8a",
"Ccr7",
"Cxcr5",
"Cxcr4",
"Ccr2",
"Pcna",
"Uhrf1"
)
split_dc <- c(
"CDP",
rep("cDC1", 4),
rep("cDC2", 13),
rep("CD", 2),
rep("Chemokines", 4),
rep("CC", 2)
)
options(repr.plot.width=2*(2.5+length(features_dc)/2.5), repr.plot.height=5)
source("bin/seurat_dea.R")
dp_1 <- dp_feature(so_m, features_dc, split_dc, group_by="varid_clusters", scale=FALSE, title="DC marker genes")
dp_2 <- dp_feature(so_m, features_dc, split_dc, group_by="varid_clusters", scale=TRUE, title="DC marker genes")
dp_1 + dp_2
`summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. `summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.
options(repr.plot.width=20, repr.plot.height=10)
varid_clusters_dc <- c("varid_clusters_m_4", "varid_clusters_m_5", "varid_clusters_m_2", "varid_clusters_m_8", "varid_clusters_m_9")
dea_dc <- dea_seurat(subset(so_m, subset=varid_clusters %in% varid_clusters_dc), ident="varid_clusters", file="result/dea/dea_dc", conserved=FALSE, only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
dea_dc_conserved <- dea_seurat(subset(so_m, subset=varid_clusters %in% varid_clusters_dc), ident="varid_clusters", file="result/dea/dea_dc_conserved", conserved=TRUE, grouping_var="treatment", only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
hm_1 <- hm_dea(dea_dc, subset(so_m, subset=varid_clusters %in% varid_clusters_dc), column_name="varid_clusters", column_order=varid_clusters_dc, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL)
hm_2 <- hm_dea(dea_dc_conserved, subset(so_m, subset=varid_clusters %in% varid_clusters_dc), column_name="varid_clusters", column_order=varid_clusters_dc, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL, conserved=TRUE)
hm_1 + hm_2
features_mo <- c(
"Ly6c2", "Csf1r", "Spi1", "Spn", # Ly6c Ly6c Csf1r Pu1 Cd43 | Canonical marker genes | Two major monocyte populations in the blood which are Ly6c lo and Ly6c high. Cd43 is low in classical and high in int and non-classical.
"Jun", "Fos", "Irf8", "Klf4", "Vcan", "Cd163", "Cd63", "S100a8", # Ccr2 | Classical monocyte | Recruitment to blood by Ccr2
"Treml4", "Cx3cr1", "Nr4a1", "Klf2", "Fcgr3", "Ifitm1", "Ifitm2", "Ifitm3", "Cdkn1c", "Mtss1", # Treml4 Cx3cr1 | Non-classical monocyte | Migrate in response to Cx3cr1
"H2-Ab1", "Cd74", "Ccr5", "Cxcl9", "Cxcl10",
"Cd4",
"Cd8a",
"Ccr7",
"Cxcr5",
"Cxcr4",
"Ccr2",
"Pcna",
"Uhrf1"
)
split_mo <- c(
rep("Canonical", 4),
rep("cMo", 8),
rep("ncMo", 10),
rep("moDC", 5),
rep("CD", 2),
rep("Chemokines", 4),
rep("CC", 2)
)
options(repr.plot.width=2*(2.5+length(features_mo)/2.5), repr.plot.height=5)
dp_1 <- dp_feature(so_m, features_mo, split_mo, group_by="varid_clusters", scale=FALSE, title="Monocyte marker")
dp_2 <- dp_feature(so_m, features_mo, split_mo, group_by="varid_clusters", scale=TRUE, title="Monocyte marker")
dp_1 + dp_2
`summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. `summarise()` has grouped output by 'varid_clusters'. You can override using the `.groups` argument. Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.
options(repr.plot.width=15, repr.plot.height=15)
varid_clusters_mo <- c("varid_clusters_m_3", "varid_clusters_m_1", "varid_clusters_m_11", "varid_clusters_m_7")
dea_mo <- dea_seurat(subset(so_m, subset=varid_clusters %in% varid_clusters_mo), ident="varid_clusters", file="result/dea/dea_mo", conserved=FALSE, only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
dea_mo_conserved <- dea_seurat(subset(so_m, subset=varid_clusters %in% varid_clusters_mo), ident="varid_clusters", file="result/dea/dea_mo_conserved", conserved=TRUE, grouping_var="treatment", only_pos=FALSE, logfc_threshold=0, min_pct=0.01, compute=FALSE)
hm_1 <- hm_dea(dea_mo, subset(so_m, subset=varid_clusters %in% varid_clusters_mo), column_name="varid_clusters", column_order=varid_clusters_mo, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL)
hm_2 <- hm_dea(dea_mo_conserved, subset(so_m, subset=varid_clusters %in% varid_clusters_mo), column_name="varid_clusters", column_order=varid_clusters_mo, row_name="ident", top=10, width=0.05, height=3, log2fc_thr=0.25, padj_thr=0.05, pct_1=0.05, pct_2=NULL, conserved=TRUE)
hm_1 + hm_2
cell_type_m <- data.frame(
ident=c(
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11
),
cell_type_fine_detail=c(
"ncMo Cd4- (1)",
"cDC2 (1)",
"cMo Ly6c lo (2)",
"cDC1 Cd8+ prolif. (1)",
"cDC1 Cd8+ (2)",
"RPM",
"ncMo Cd4+ (2)",
"cDC2 Ccr7+ (3)",
"cDC2 prolif. (2)",
"PreRPM",
"cMo Ly6c hi (1)"
),
cell_type_fine=c(
"ncMo (1)",
"cDC2 (1)",
"cMo (2)",
"cDC1 (1)",
"cDC1 (2)",
"RPM",
"ncMo (2)",
"cDC2 (3)",
"cDC2 (2)",
"PreRPM",
"cMo (1)"
),
cell_type_main=c(
"ncMo",
"cDC2",
"cMo",
"cDC1",
"cDC1",
"RPM",
"ncMo",
"cDC2",
"cDC2",
"PreRPM",
"cMo"
)
)
cell_type_m <- dplyr::mutate(cell_type_m, ident=paste0("varid_clusters_m_", ident)) %>%
dplyr::left_join(dplyr::select(so_m@meta.data, seurat_clusters, varid_clusters, cell_id), ., by=c("varid_clusters"="ident")) %>%
column_to_rownames("cell_id") %>%
dplyr::mutate(ident=varid_clusters)
cell_type <- rbind(cell_type_eb, cell_type_mpp, cell_type_m) %>% dplyr::select(cell_type_main, cell_type_fine, cell_type_fine_detail, varid_clusters, ident)
so <- AddMetaData(so, cell_type)
unique(dplyr::select(cell_type, cell_type_main, cell_type_fine, cell_type_fine_detail, ident)) %>% write.csv("test.csv")
options(repr.plot.width=10, repr.plot.height=5)
# Source files
source("plotting_global.R")
source("bin/seurat_qc.R")
source("bin/cell_type.R")
dplot_1 <- dplot(so, reduction="umap", group_by="cell_type_fine_detail", alpha=1, pt_size=0.5) +
scale_color_manual(values=color[["cell_type_fine_detail"]]) +
theme(
legend.title=element_blank(),
plot.title=element_blank()
)
dplot_2 <- dplot(so, reduction="umap", group_by="cell_type_main", alpha=1, pt_size=0.5) +
scale_color_manual(values=color[["cell_type_main"]]) +
theme(
legend.title=element_blank(),
plot.title=element_blank()
)
dplot_1 + dplot_2
saveRDS(so, so_pp_file)
write.csv(so@meta.data, "data/object/components/meta.csv")
write.csv(rownames(so), "data/object/genes.csv", row.names=FALSE)
write.csv(colnames(so), "data/object/cells.csv", row.names=FALSE)
write.csv(so@reductions$umap@cell.embeddings, "data/object/components/umap.csv")
# Store data as h5ad
adata <- import("anndata", as="adata", convert=FALSE)
pd <- import("pandas", as="pd", convert=FALSE)
np <- import("numpy", as="np", convert=FALSE)
# Transform dgCMatrix to
X <- GetAssayData(so, assay="RNA", slot="counts") %>% as.matrix() %>% t()
X <- np$array(X, dtype=np$int32)
adata <- adata$AnnData(X=X, obs=so@meta.data)
adata$var_names <- rownames(GetAssayData(so, assay="RNA", slot="counts"))
adata$raw <- adata
adata$write_h5ad(h5ad_pp_file)
None
sessionInfo()
R version 4.1.0 (2021-05-18) Platform: x86_64-conda-linux-gnu (64-bit) Running under: Red Hat Enterprise Linux 8.5 (Ootpa) Matrix products: default BLAS/LAPACK: /home/fdeckert/bin/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.15.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] grid stats4 stats graphics grDevices utils datasets [8] methods base other attached packages: [1] patchwork_1.1.1 RColorBrewer_1.1-3 [3] reticulate_1.24 ggplotify_0.1.0 [5] viridis_0.6.2 viridisLite_0.4.0 [7] circlize_0.4.14 ComplexHeatmap_2.10.0 [9] openxlsx_4.2.5 forcats_0.5.1 [11] stringr_1.4.0 dplyr_1.0.9 [13] purrr_0.3.4 readr_2.1.2 [15] tidyr_1.2.0 tibble_3.1.7 [17] ggplot2_3.3.6 tidyverse_1.3.1 [19] fgsea_1.20.0 msigdbr_7.5.1 [21] clusterExperiment_2.14.0 tradeSeq_1.8.0 [23] SingleR_1.8.1 RaceID_0.2.6 [25] Matrix_1.4-1 monocle3_1.0.0 [27] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0 [29] GenomicRanges_1.46.1 GenomeInfoDb_1.30.1 [31] IRanges_2.28.0 S4Vectors_0.32.4 [33] MatrixGenerics_1.6.0 matrixStats_0.62.0 [35] Biobase_2.54.0 BiocGenerics_0.40.0 [37] SeuratWrappers_0.3.0 sp_1.4-7 [39] SeuratObject_4.1.0 Seurat_4.1.1 loaded via a namespace (and not attached): [1] rsvd_1.0.5 svglite_2.1.0 [3] ica_1.0-2 zinbwave_1.16.0 [5] foreach_1.5.2 lmtest_0.9-40 [7] crayon_1.5.1 spatstat.core_2.4-2 [9] MASS_7.3-57 rhdf5filters_1.6.0 [11] nlme_3.1-157 backports_1.4.1 [13] reprex_2.0.1 rlang_1.0.6 [15] XVector_0.34.0 ROCR_1.0-11 [17] readxl_1.4.0 irlba_2.3.5 [19] limma_3.50.3 phylobase_0.8.10 [21] runner_0.4.2 BiocParallel_1.28.3 [23] rjson_0.2.21 snowfall_1.84-6.1 [25] bit64_4.0.5 glue_1.6.2 [27] harmony_0.1.0 pheatmap_1.0.12 [29] rngtools_1.5.2 sctransform_0.3.3 [31] parallel_4.1.0 spatstat.sparse_3.0-0 [33] AnnotationDbi_1.56.2 coop_0.6-3 [35] spatstat.geom_3.0-5 haven_2.5.0 [37] tidyselect_1.1.2 fitdistrplus_1.1-8 [39] XML_3.99-0.9 zoo_1.8-10 [41] xtable_1.8-4 magrittr_2.0.3 [43] evaluate_0.15 cli_3.4.1 [45] zlibbioc_1.40.0 rstudioapi_0.13 [47] miniUI_0.1.1.1 rpart_4.1.16 [49] fastmatch_1.1-3 locfdr_1.1-8 [51] shiny_1.7.1 xfun_0.30 [53] BiocSingular_1.10.0 clue_0.3-60 [55] askpass_1.1 cluster_2.1.3 [57] pbdZMQ_0.3-7 KEGGREST_1.34.0 [59] ggrepel_0.9.1 ape_5.6-2 [61] listenv_0.8.0 Biostrings_2.62.0 [63] png_0.1-7 permute_0.9-7 [65] future_1.25.0 withr_2.5.0 [67] bitops_1.0-7 plyr_1.8.7 [69] cellranger_1.1.0 pillar_1.8.1 [71] GlobalOptions_0.1.2 cachem_1.0.6 [73] fs_1.5.2 kernlab_0.9-30 [75] scatterplot3d_0.3-41 GetoptLong_1.0.5 [77] DelayedMatrixStats_1.16.0 vctrs_0.5.1 [79] ellipsis_0.3.2 generics_0.1.2 [81] NMF_0.24.0 tools_4.1.0 [83] rncl_0.8.6 munsell_0.5.0 [85] DelayedArray_0.20.0 fastmap_1.1.0 [87] compiler_4.1.0 abind_1.4-5 [89] httpuv_1.6.5 pkgmaker_0.32.2 [91] plotly_4.10.0 rgeos_0.5-9 [93] GenomeInfoDbData_1.2.7 gridExtra_2.3 [95] edgeR_3.36.0 lattice_0.20-45 [97] deldir_1.0-6 utf8_1.2.2 [99] later_1.3.0 jsonlite_1.8.0 [101] scales_1.2.0 princurve_2.1.6 [103] ScaledMatrix_1.2.0 pbapply_1.5-0 [105] sparseMatrixStats_1.6.0 genefilter_1.76.0 [107] lazyeval_0.2.2 promises_1.2.0.1 [109] doParallel_1.0.17 R.utils_2.11.0 [111] goftest_1.2-3 spatstat.utils_3.0-1 [113] rmarkdown_2.14 cowplot_1.1.1 [115] webshot_0.5.3 Rtsne_0.16 [117] softImpute_1.4-1 uwot_0.1.11 [119] igraph_1.3.1 HDF5Array_1.22.1 [121] survival_3.3-1 systemfonts_1.0.4 [123] htmltools_0.5.2 memoise_2.0.1 [125] locfit_1.5-9.5 quadprog_1.5-8 [127] digest_0.6.29 assertthat_0.2.1 [129] mime_0.12 repr_1.1.4 [131] registry_0.5-1 RSQLite_2.2.13 [133] yulab.utils_0.0.4 future.apply_1.9.0 [135] som_0.3-5.1 remotes_2.4.2 [137] data.table_1.14.2 blob_1.2.3 [139] vegan_2.6-2 R.oo_1.24.0 [141] RNeXML_2.4.7 labeling_0.4.2 [143] splines_4.1.0 Rhdf5lib_1.16.0 [145] Cairo_1.6-0 RCurl_1.98-1.6 [147] broom_0.8.0 hms_1.1.1 [149] modelr_0.1.8 rhdf5_2.38.1 [151] colorspace_2.0-3 base64enc_0.1-3 [153] BiocManager_1.30.17 shape_1.4.6 [155] Rcpp_1.0.8.3 RANN_2.6.1 [157] fansi_1.0.3 tzdb_0.3.0 [159] parallelly_1.31.1 IRdisplay_1.1 [161] R6_2.5.1 ggridges_0.5.3 [163] lifecycle_1.0.3 zip_2.2.0 [165] leiden_0.3.10 howmany_0.3-1 [167] RcppAnnoy_0.0.19 iterators_1.0.14 [169] htmlwidgets_1.5.4 umap_0.2.8.0 [171] beachmat_2.10.0 polyclip_1.10-0 [173] propr_4.2.6 gridGraphics_0.5-1 [175] rvest_1.0.2 mgcv_1.8-40 [177] globals_0.14.0 openssl_2.0.0 [179] spatstat.random_3.1-3 lle_1.1 [181] slingshot_2.2.1 progressr_0.10.0 [183] codetools_0.2-18 lubridate_1.8.0 [185] FNN_1.1.3 randomForest_4.7-1 [187] prettyunits_1.1.1 dbplyr_2.1.1 [189] gridBase_0.4-7 RSpectra_0.16-1 [191] R.methodsS3_1.8.1 gtable_0.3.0 [193] DBI_1.1.2 highr_0.9 [195] tensor_1.5 httr_1.4.3 [197] KernSmooth_2.23-20 stringi_1.7.6 [199] progress_1.2.2 farver_2.1.0 [201] reshape2_1.4.4 uuid_1.1-0 [203] annotate_1.72.0 magick_2.7.3 [205] xml2_1.3.3 IRkernel_1.3 [207] kableExtra_1.3.4 BiocNeighbors_1.12.0 [209] ade4_1.7-19 scattermore_0.8 [211] FateID_0.2.1 bit_4.0.4 [213] spatstat.data_3.0-0 pkgconfig_2.0.3 [215] babelgene_22.3 TrajectoryUtils_1.2.0 [217] knitr_1.39